#include "spglib.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static void test_spg_get_symmetry(void);
static void test_spg_get_symmetry_with_collinear_spin(void);
static void test_spg_get_multiplicity(void);
static void test_spg_find_primitive_BCC(void);
static void test_spg_find_primitive_corundum(void);
static void test_spg_standardize_cell_BCC(void);
static void test_spg_standardize_cell_BCC_prim(void);
static void test_spg_standardize_cell_corundum(void);
static int sub_spg_standardize_cell(double lattice[3][3],
				    double position[][3],
				    int types[],
				    const int num_atom,
				    const double symprec,
				    const int to_primitive,
				    const int no_idealize);
static void test_spg_get_international(void);
static void test_spg_get_schoenflies(void);
static void test_spg_get_spacegroup_type(void);
static void test_spg_get_symmetry_from_database(void);
static void test_spg_refine_cell_BCC(void);
static void test_spg_get_dataset(void);
static void test_spg_get_ir_reciprocal_mesh(void);
static void test_spg_get_stabilized_reciprocal_mesh(void);
static void test_spg_relocate_BZ_grid_address(void);
static void test_spg_get_error_message(void);
static void show_spg_dataset(double lattice[3][3],
			     const double origin_shift[3],
			     double position[][3],
			     const int num_atom,
			     const int types[]);
static void show_cell(double lattice[3][3],
		      double position[][3],
		      const int types[],
		      const int num_atom);

int main(void)
{
  test_spg_find_primitive_BCC();
  test_spg_find_primitive_corundum();
  test_spg_standardize_cell_BCC();
  test_spg_standardize_cell_BCC_prim();
  test_spg_standardize_cell_corundum();
  test_spg_get_multiplicity();
  test_spg_get_symmetry();
  test_spg_get_symmetry_with_collinear_spin();
  test_spg_get_international();
  test_spg_get_schoenflies();
  test_spg_get_spacegroup_type();
  test_spg_get_symmetry_from_database();
  test_spg_refine_cell_BCC();
  test_spg_get_dataset();
  test_spg_get_ir_reciprocal_mesh();
  test_spg_get_stabilized_reciprocal_mesh();
  test_spg_relocate_BZ_grid_address();
  test_spg_get_error_message();

  return 0;
}

static void test_spg_find_primitive_BCC(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 4}};
  double position[][3] = {
    {0, 0, 0},
    {0.5, 0.5, 0.5}
  };
  int types[] = {1, 1};
  int i, num_atom = 2, num_primitive_atom;
  double symprec = 1e-5;
  
  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_find_primitive (BCC unitcell --> primitive) ***:\n");
  num_primitive_atom = spg_find_primitive(lattice, position, types, num_atom, symprec);
  if (num_primitive_atom == 0) {
    printf("Primitive cell was not found.\n");
  } else { 
    show_cell(lattice, position, types, num_primitive_atom);
  }
}

static void test_spg_find_primitive_corundum(void)
{
  double lattice[3][3] = {{4.8076344022756095, -2.4038172011378047, 0},
			  {0, 4.1635335244786962, 0},
			  {0, 0, 13.1172699198127543}};
  double position[][3] = {
    {0.0000000000000000, 0.0000000000000000, 0.3521850942289043},
    {0.6666666666666643, 0.3333333333333357, 0.6855184275622400},
    {0.3333333333333357, 0.6666666666666643, 0.0188517608955686},
    {0.0000000000000000, 0.0000000000000000, 0.6478149057711028},
    {0.6666666666666643, 0.3333333333333357, 0.9811482391044314},
    {0.3333333333333357, 0.6666666666666643, 0.3144815724377600},
    {0.0000000000000000, 0.0000000000000000, 0.1478149057710957},
    {0.6666666666666643, 0.3333333333333357, 0.4811482391044314},
    {0.3333333333333357, 0.6666666666666643, 0.8144815724377600},
    {0.0000000000000000, 0.0000000000000000, 0.8521850942288972},
    {0.6666666666666643, 0.3333333333333357, 0.1855184275622400},
    {0.3333333333333357, 0.6666666666666643, 0.5188517608955686},
    {0.3061673906454899, 0.0000000000000000, 0.2500000000000000},
    {0.9728340573121541, 0.3333333333333357, 0.5833333333333357},
    {0.6395007239788255, 0.6666666666666643, 0.9166666666666643},
    {0.6938326093545102, 0.0000000000000000, 0.7500000000000000},
    {0.3604992760211744, 0.3333333333333357, 0.0833333333333357},
    {0.0271659426878458, 0.6666666666666643, 0.4166666666666643},
    {0.0000000000000000, 0.3061673906454899, 0.2500000000000000},
    {0.6666666666666643, 0.6395007239788255, 0.5833333333333357},
    {0.3333333333333357, 0.9728340573121541, 0.9166666666666643},
    {0.0000000000000000, 0.6938326093545102, 0.7500000000000000},
    {0.6666666666666643, 0.0271659426878458, 0.0833333333333357},
    {0.3333333333333357, 0.3604992760211744, 0.4166666666666643},
    {0.6938326093545102, 0.6938326093545102, 0.2500000000000000},
    {0.3604992760211744, 0.0271659426878458, 0.5833333333333357},
    {0.0271659426878458, 0.3604992760211744, 0.9166666666666643},
    {0.3061673906454899, 0.3061673906454899, 0.7500000000000000},
    {0.9728340573121541, 0.6395007239788255, 0.0833333333333357},
    {0.6395007239788255, 0.9728340573121541, 0.4166666666666643},
  };
  int types[30];
  int i, num_atom = 30, num_primitive_atom;
  double symprec = 1e-5;
  
  for (i = 0; i < 12; i++) {
    types[i] = 1;
  }
  for (i = 12; i < 30; i++) {
    types[i] = 2;
  }

  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_find_primitive (Corundum) ***:\n");
  num_primitive_atom = spg_find_primitive(lattice, position, types, num_atom, symprec);
  if (num_primitive_atom == 0) {
    printf("Primitive cell was not found.\n");
  } else { 
    show_cell(lattice, position, types, num_primitive_atom);
  }
}

static void test_spg_refine_cell_BCC(void)
{
  double lattice[3][3] = {{0, 2, 2}, {2, 0, 2}, {2, 2, 0}};

  /* 4 times larger memory space must be prepared. */
  double position[4][3];
  int types[4];
  
  int i, num_atom_bravais, num_atom = 1;
  double symprec = 1e-5;

  position[0][0] = 0;
  position[0][1] = 0;
  position[0][2] = 0;
  types[0] = 1;
  
  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_refine_cell ***:\n");
  num_atom_bravais = spg_refine_cell( lattice,
				      position,
				      types,
				      num_atom,
				      symprec );
  show_cell(lattice, position, types, num_atom_bravais);
}

static void test_spg_standardize_cell_BCC(void)
{
  double lattice[3][3] = {{3.97, 0, 0}, {0, 4.03, 0}, {0, 0, 4.0}};
  double position[][3] = {
    {0.002, 0, 0},
    {0.5, 0.5001, 0.5}
  };
  int types[] = {1, 1};
  int i, j, k, num_atom = 2, num_primitive_atom;
  double symprec = 1e-1;
  
  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_standardize_cell (BCC unitcell) ***:\n");
  printf("------------------------------------------------------\n");
  for (j = 0; j < 2; j++) {
    for (k = 0; k < 2; k++) {
      sub_spg_standardize_cell(lattice,
			       position,
			       types,
			       num_atom,
			       symprec,
			       j,
			       k);
      printf("------------------------------------------------------\n");
    }
  }
}

static void test_spg_standardize_cell_BCC_prim(void)
{
  double lattice[3][3] = {{-2.01, 2, 2}, {2, -2.02, 2}, {2, 2, -2.03}};
  double position[][3] = {
    {0.002, 0, 0},
  };
  int types[] = {1};
  int i, j, k, num_atom = 1, num_primitive_atom;
  double symprec = 1e-1;
  
  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_standardize_cell (BCC primitive) ***:\n");
  printf("------------------------------------------------------\n");
  for (j = 0; j < 2; j++) {
    for (k = 0; k < 2; k++) {
      sub_spg_standardize_cell(lattice,
			       position,
			       types,
			       num_atom,
			       symprec,
			       j,
			       k);
      printf("------------------------------------------------------\n");
    }
  }
}

static void test_spg_standardize_cell_corundum(void)
{
  double lattice[3][3] = {{4.8076344022756095, -2.4038172011378047, 0},
			  {0, 4.1635335244786962, 0},
			  {0, 0, 13.1172699198127543}};
  double position[][3] = {
    {0.0000000000000000, 0.0000000000000000, 0.3521850942289043},
    {0.6666666666666643, 0.3333333333333357, 0.6855184275622400},
    {0.3333333333333357, 0.6666666666666643, 0.0188517608955686},
    {0.0000000000000000, 0.0000000000000000, 0.6478149057711028},
    {0.6666666666666643, 0.3333333333333357, 0.9811482391044314},
    {0.3333333333333357, 0.6666666666666643, 0.3144815724377600},
    {0.0000000000000000, 0.0000000000000000, 0.1478149057710957},
    {0.6666666666666643, 0.3333333333333357, 0.4811482391044314},
    {0.3333333333333357, 0.6666666666666643, 0.8144815724377600},
    {0.0000000000000000, 0.0000000000000000, 0.8521850942288972},
    {0.6666666666666643, 0.3333333333333357, 0.1855184275622400},
    {0.3333333333333357, 0.6666666666666643, 0.5188517608955686},
    {0.3061673906454899, 0.0000000000000000, 0.2500000000000000},
    {0.9728340573121541, 0.3333333333333357, 0.5833333333333357},
    {0.6395007239788255, 0.6666666666666643, 0.9166666666666643},
    {0.6938326093545102, 0.0000000000000000, 0.7500000000000000},
    {0.3604992760211744, 0.3333333333333357, 0.0833333333333357},
    {0.0271659426878458, 0.6666666666666643, 0.4166666666666643},
    {0.0000000000000000, 0.3061673906454899, 0.2500000000000000},
    {0.6666666666666643, 0.6395007239788255, 0.5833333333333357},
    {0.3333333333333357, 0.9728340573121541, 0.9166666666666643},
    {0.0000000000000000, 0.6938326093545102, 0.7500000000000000},
    {0.6666666666666643, 0.0271659426878458, 0.0833333333333357},
    {0.3333333333333357, 0.3604992760211744, 0.4166666666666643},
    {0.6938326093545102, 0.6938326093545102, 0.2500000000000000},
    {0.3604992760211744, 0.0271659426878458, 0.5833333333333357},
    {0.0271659426878458, 0.3604992760211744, 0.9166666666666643},
    {0.3061673906454899, 0.3061673906454899, 0.7500000000000000},
    {0.9728340573121541, 0.6395007239788255, 0.0833333333333357},
    {0.6395007239788255, 0.9728340573121541, 0.4166666666666643},
  };
  int types[30];
  int i, j, k, num_atom = 30;
  double symprec = 1e-5;
  
  for (i = 0; i < 12; i++) {
    types[i] = 1;
  }
  for (i = 12; i < 30; i++) {
    types[i] = 2;
  }
  
  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_standardize_cell (Corundum) ***:\n");
  printf("------------------------------------------------------\n");
  for (j = 0; j < 2; j++) {
    for (k = 0; k < 2; k++) {
      sub_spg_standardize_cell(lattice,
			       position,
			       types,
			       num_atom,
			       symprec,
			       j,
			       k);
      printf("------------------------------------------------------\n");
    }
  }
}

static int sub_spg_standardize_cell(double lattice[3][3],
				    double position[][3],
				    int types[],
				    const int num_atom,
				    const double symprec,
				    const int to_primitive,
				    const int no_idealize)
{
  int i, num_primitive_atom;
  double lat[3][3], pos[num_atom][3];
  int typ[num_atom];

  for (i = 0; i < 3; i++) {
    lat[i][0] = lattice[i][0];
    lat[i][1] = lattice[i][1];
    lat[i][2] = lattice[i][2];
  }

  for (i = 0; i < num_atom; i++) {
    pos[i][0] = position[i][0];
    pos[i][1] = position[i][1];
    pos[i][2] = position[i][2];
    typ[i] = types[i];
  }
  
  /* lattice, position, and types are overwirtten. */
  num_primitive_atom = spg_standardize_cell(lat,
					    pos,
					    typ,
					    num_atom,
					    to_primitive,
					    no_idealize,
					    symprec);
  printf("VASP POSCAR format: ");
  if (to_primitive == 0) {
    printf("to_primitive=0 and ");
  } else {
    printf("to_primitive=1 and ");
  }

  if (no_idealize == 0) {
    printf("no_idealize=0\n");
  } else {
    printf("no_idealize=1\n");
  }
  printf("1.0\n");
  for (i = 0; i < 3; i++) {
    printf("%f %f %f\n", lat[0][i], lat[1][i], lat[2][i]);
  }
  printf("%d\n", num_primitive_atom);
  printf("Direct\n");
  for (i = 0; i < num_primitive_atom; i++) {
    printf("%f %f %f\n", pos[i][0], pos[i][1], pos[i][2]);
  }
}

static void test_spg_get_international(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.5},
      {0.3, 0.3, 0},
      {0.7, 0.7, 0},
      {0.2, 0.8, 0.5},
      {0.8, 0.2, 0.5},
    };
  int types[] = {1, 1, 2, 2, 2, 2};
  int num_spg, num_atom = 6;
  char symbol[21];
  
  num_spg = spg_get_international(symbol, lattice, position, types, num_atom, 1e-5);
  printf("*** Example of spg_get_international ***:\n");
  if ( num_spg > 0 ) {
    printf("%s (%d)\n", symbol, num_spg);
  } else {
    printf("Space group could not be found.\n");
  }
}

static void test_spg_get_schoenflies(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.5},
      {0.3, 0.3, 0},
      {0.7, 0.7, 0},
      {0.2, 0.8, 0.5},
      {0.8, 0.2, 0.5},
    };
  int types[] = {1, 1, 2, 2, 2, 2};
  int num_atom = 6;
  char symbol[7];
  
  spg_get_schoenflies(symbol, lattice, position, types, num_atom, 1e-5);
  printf("*** Example of spg_get_schoenflies ***:\n");
  printf("Schoenflies: %s\n", symbol);
}

static void test_spg_get_spacegroup_type(void)
{
  SpglibSpacegroupType spgtype;
  spgtype = spg_get_spacegroup_type(446);

  printf("*** Example of spg_get_spacegroup_type ***:\n");
  printf("Number:        %d\n", spgtype.number);
  printf("Schoenflies:   %s\n", spgtype.schoenflies);
  printf("International: %s\n", spgtype.international);
  printf("International: %s\n", spgtype.international_full);
  printf("International: %s\n", spgtype.international_short);
  printf("Hall symbol:   %s\n", spgtype.hall_symbol);
  
}

static void test_spg_get_symmetry_from_database(void)
{
  int rotations[192][3][3];
  double translations[192][3];
  int max_size, i, j, size;
  
  size = spg_get_symmetry_from_database(rotations,
					translations,
					460);
  
  printf("*** Example of spg_get_symmetry_from_database ***:\n");
  for (i = 0; i < size; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++) {
      printf("%2d %2d %2d\n",
	     rotations[i][j][0], rotations[i][j][1], rotations[i][j][2]);
    }
    printf("%f %f %f\n",
	   translations[i][0], translations[i][1], translations[i][2]);
  }
}

static void test_spg_get_multiplicity(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 4}};
  double position[][3] = {
    {0, 0, 0},
    {0.5, 0.5, 0.5}
  };
  int types[] = {1, 2};
  int num_atom = 2;
  int size;

  size = spg_get_multiplicity(lattice, position, types, num_atom, 1e-5);

  printf("*** Example of spg_get_multiplicity ***:\n");
  printf("Number of symmetry operations: %d\n", size);
}


static void test_spg_get_symmetry(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.25},
      {0.3, 0.3, 0},
      {0.7, 0.7, 0},
      {0.2, 0.8, 0.25},
      {0.8, 0.2, 0.25},
      {0, 0, 0.5},
      {0.5, 0.5, 0.75},
      {0.3, 0.3, 0.5},
      {0.7, 0.7, 0.5},
      {0.2, 0.8, 0.75},
      {0.8, 0.2, 0.75}
    };
  int types[] = {1, 1, 2, 2, 2, 2, 1, 1, 2, 2, 2, 2};
  int num_atom = 12;
  int max_size = 50;
  int i, j, size;
  int rotation[max_size][3][3];
  double translation[max_size][3];

  double origin_shift[3] = {0.1, 0.1, 0};
  for (i = 0; i < num_atom; i++) {
    for (j = 0; j < 3; j++) {
      position[i][j] += origin_shift[j];
    }
  }

  printf("*** Example of spg_get_symmetry (Rutile two unit cells) ***:\n");
  size = spg_get_symmetry(rotation,
			  translation,
			  max_size,
			  lattice,
			  position,
			  types,
			  num_atom,
			  1e-5);
  for (i = 0; i < size; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++)
      printf("%2d %2d %2d\n", rotation[i][j][0], rotation[i][j][1], rotation[i][j][2]);
    printf("%f %f %f\n",
	   translation[i][0], translation[i][1], translation[i][2]);
  }

}

static void test_spg_get_symmetry_with_collinear_spin(void) {
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 4}};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.5}
    };	
  int types[] = {1, 1};
  int equivalent_atoms[2];
  double spins[2];
  int num_atom = 2;
  int max_size = 300;
  int i, j, size;
  int rotation[max_size][3][3];
  double translation[max_size][3];

  printf("*** Example of spg_get_symmetry_with_spin (BCC ferro) ***:\n");
  spins[0] = 1;
  spins[1] = 1;
  size = spg_get_symmetry_with_collinear_spin(rotation,
					      translation,
					      equivalent_atoms,
					      max_size,
					      lattice,
					      position,
					      types,
					      spins,
					      num_atom,
					      1e-5);
  for (i = 0; i < size; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++) {
      printf("%2d %2d %2d\n",
	     rotation[i][j][0], rotation[i][j][1], rotation[i][j][2]);
    }
    printf("%f %f %f\n", translation[i][0], translation[i][1],
           translation[i][2]);
  }

  printf("*** Example of spg_get_symmetry_with_spin (BCC antiferro) ***:\n");
  spins[0] = 1;
  spins[1] = -1;
  size = spg_get_symmetry_with_collinear_spin(rotation,
					      translation,
					      equivalent_atoms,
					      max_size,
					      lattice,
					      position,
					      types,
					      spins,
					      num_atom,
					      1e-5);
  for (i = 0; i < size; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++) {
      printf("%2d %2d %2d\n",
	     rotation[i][j][0], rotation[i][j][1], rotation[i][j][2]);
    }
    printf("%f %f %f\n", translation[i][0], translation[i][1],
           translation[i][2]);
  }

  printf("*** Example of spg_get_symmetry_with_spin (BCC broken spin) ***:\n");
  spins[0] = 1;
  spins[1] = 2;
  size = spg_get_symmetry_with_collinear_spin(rotation,
					      translation,
					      equivalent_atoms,
					      max_size,
					      lattice,
					      position,
					      types,
					      spins,
					      num_atom,
					      1e-5);
  for (i = 0; i < size; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++) {
      printf("%2d %2d %2d\n",
	     rotation[i][j][0], rotation[i][j][1], rotation[i][j][2]);
    }
    printf("%f %f %f\n", translation[i][0], translation[i][1],
           translation[i][2]);
  }
}


static void test_spg_get_dataset(void)
{
  printf("*** Example of spg_get_dataset (Rutile two unit cells) ***:\n");
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
  double origin_shift[3] = {0.1, 0.1, 0};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.25},
      {0.3, 0.3, 0},
      {0.7, 0.7, 0},
      {0.2, 0.8, 0.25},
      {0.8, 0.2, 0.25},
      {0, 0, 0.5},
      {0.5, 0.5, 0.75},
      {0.3, 0.3, 0.5},
      {0.7, 0.7, 0.5},
      {0.2, 0.8, 0.75},
      {0.8, 0.2, 0.75}
    };
  int types[] = {1, 1, 2, 2, 2, 2, 1, 1, 2, 2, 2, 2};
  int num_atom = 12;

  show_spg_dataset(lattice, origin_shift, position, num_atom, types);

  double lattice_2[3][3] = {{3.7332982433264039, -1.8666491216632011, 0},
			    {0, 3.2331311186244847, 0},
			    {0, 0, 6.0979971306362799}};
  double origin_shift_2[3] = {0.1, 0.1, 0};
  double position_2[][3] =
    {
      {0, 0, 0},
      {1.0 / 3, 2.0 / 3, 0.4126},
      {1.0 / 3, 2.0 / 3, 0.776},
      {2.0 / 3, 1.0 / 3, 0.2542},
    };
  int types_2[] = {1, 2, 3, 3};
  int num_atom_2 = 4;

  show_spg_dataset(lattice_2, origin_shift_2, position_2, num_atom_2, types_2);
}

static void test_spg_get_ir_reciprocal_mesh(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.5},
      {0.3, 0.3, 0},
      {0.7, 0.7, 0},
      {0.2, 0.8, 0.5},
      {0.8, 0.2, 0.5},
    };
  int types[] = {1, 1, 2, 2, 2, 2};
  int num_atom = 6;
  int m = 40;
  int mesh[] = {m, m, m};
  int is_shift[] = {1, 1, 1};
  int grid_address[m * m * m][3];
  int grid_mapping_table[m * m * m];

  printf("*** Example of spg_get_ir_reciprocal_mesh of Rutile structure ***:\n");

  int num_ir = spg_get_ir_reciprocal_mesh(grid_address,
					  grid_mapping_table,
					  mesh,
					  is_shift,
					  1,
					  lattice,
					  position,
					  types,
					  num_atom,
					  1e-5);

  printf("Number of irreducible k-points of Rutile with\n");
  printf("40x40x40 Monkhorst-Pack mesh is %d (4200).\n", num_ir);
}

static void test_spg_get_stabilized_reciprocal_mesh(void)
{
  SpglibDataset *dataset;
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
  double position[][3] =
    {
      {0, 0, 0},
      {0.5, 0.5, 0.5},
      {0.3, 0.3, 0},
      {0.7, 0.7, 0},
      {0.2, 0.8, 0.5},
      {0.8, 0.2, 0.5},
    };
  int types[] = {1, 1, 2, 2, 2, 2};
  int num_atom = 6;

  dataset = spg_get_dataset(lattice,
			    position,
			    types,
			    num_atom,
			    1e-5);

  
  int m = 40;
  int mesh[] = {m, m, m};
  int is_shift[] = {1, 1, 1};
  int grid_address[m * m * m][3];
  int grid_mapping_table[m * m * m];
  double q[] = {0, 0.5, 0.5};

  printf("*** Example of spg_get_stabilized_reciprocal_mesh of Rutile structure ***:\n");

  int num_ir = spg_get_stabilized_reciprocal_mesh(grid_address,
						  grid_mapping_table,
						  mesh,
						  is_shift,
						  1,
						  dataset->n_operations,
						  dataset->rotations,
						  1,
						  (double(*)[3])q);
  spg_free_dataset(dataset);

  printf("Number of irreducible k-points stabilized by q=(0, 1/2, 1/2) of Rutile with\n");
  printf("40x40x40 Monkhorst-Pack mesh is %d (8000).\n", num_ir);
}

static void test_spg_relocate_BZ_grid_address(void)
{
  double rec_lattice[3][3] = {{-0.17573761,  0.17573761,  0.17573761},
			      { 0.17573761, -0.17573761,  0.17573761},
			      { 0.17573761,  0.17573761, -0.17573761}};
  int rotations[][3][3] = {{{1, 0, 0},
			    {0, 1, 0},
			    {0, 0, 1}}};

  
  int m = 40;
  int mesh[] = {m, m, m};
  int is_shift[] = {0, 0, 0};
  int bz_grid_address[(m + 1) * (m + 1) * (m + 1)][3];
  int bz_map[m * m * m * 8];
  int grid_address[m * m * m][3];
  int grid_mapping_table[m * m * m];
  double q[] = {0, 0, 0};

  int num_ir = spg_get_stabilized_reciprocal_mesh(grid_address,
						  grid_mapping_table,
						  mesh,
						  is_shift,
						  1,
						  1,
						  rotations,
						  1,
						  (double(*)[3])q);
  

  printf("*** Example of spg_relocate_BZ_grid_address of NaCl structure ***:\n");

  int num_q = spg_relocate_BZ_grid_address(bz_grid_address,
					   bz_map,
					   grid_address,
					   mesh,
					   rec_lattice,
					   is_shift);

  printf("Number of k-points of NaCl Brillouin zone\n");
  printf("with Gamma-centered 40x40x40 Monkhorst-Pack mesh is %d (65861).\n", num_q);
}

static void show_spg_dataset(double lattice[3][3],
			     const double origin_shift[3],
			     double position[][3],
			     const int num_atom,
			     const int types[])
{
  SpglibDataset *dataset;
  char ptsymbol[6];
  int pt_trans_mat[3][3];

  int i, j, size;
  const char *wl = "abcdefghijklmnopqrstuvwxyz";

  for ( i = 0; i < num_atom; i++ ) {
    for ( j = 0; j < 3; j++ ) {
      position[i][j] += origin_shift[j];
    }
  }

  dataset = spg_get_dataset(lattice,
			    position,
			    types,
			    num_atom,
			    1e-5);
  
  printf("International: %s (%d)\n", dataset->international_symbol, dataset->spacegroup_number );
  printf("Hall symbol:   %s\n", dataset->hall_symbol );
  spg_get_pointgroup(ptsymbol,
		     pt_trans_mat,
		     dataset->rotations,
		     dataset->n_operations);
  printf("Point group:   %s\n", ptsymbol);
  printf("Transformation matrix:\n");
  for ( i = 0; i < 3; i++ ) {
    printf("%f %f %f\n",
	   dataset->transformation_matrix[i][0],
	   dataset->transformation_matrix[i][1],
	   dataset->transformation_matrix[i][2]);
  }
  printf("Wyckoff letters:\n");
  for ( i = 0; i < dataset->n_atoms; i++ ) {
    printf("%c ", wl[dataset->wyckoffs[i]]);
  }
  printf("\n");
  printf("Equivalent atoms:\n");
  for (i = 0; i < dataset->n_atoms; i++) {
    printf("%d ", dataset->equivalent_atoms[i]);
  }
  printf("\n");
  
  for (i = 0; i < dataset->n_operations; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++) {
      printf("%2d %2d %2d\n",
	     dataset->rotations[i][j][0],
	     dataset->rotations[i][j][1],
	     dataset->rotations[i][j][2]);
    }
    printf("%f %f %f\n",
	   dataset->translations[i][0],
	   dataset->translations[i][1],
	   dataset->translations[i][2]);
  }

  spg_free_dataset(dataset);

}

static void show_cell(double lattice[3][3],
		      double position[][3],
		      const int types[],
		      const int num_atom)
{
  int i;

  printf("Lattice parameter:\n");
  for (i = 0; i < 3; i++) {
    printf("%f %f %f\n", lattice[0][i], lattice[1][i], lattice[2][i]);
  }
  printf("Atomic positions:\n");
  for (i = 0; i < num_atom; i++) {
    printf("%d: %f %f %f\n",
	   types[i], position[i][0], position[i][1], position[i][2]);
  }
}

static void test_spg_get_error_message(void)
{
  double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 4}};
  double position[][3] = {
    {0, 0, 0},
    {0.5, 0.5, 0.5},
    {0.5, 0.5, 0.5}
  };
  int types[] = {1, 1, 1};
  int i, num_atom = 3, num_primitive_atom;
  double symprec = 1e-5;
  SpglibError error;
  
  
  /* lattice, position, and types are overwirtten. */
  printf("*** Example of spg_get_error_message ***:\n");
  num_primitive_atom = spg_find_primitive(lattice, position, types, num_atom, symprec);
  if (num_primitive_atom == 0) {
    printf("Primitive cell was not found.\n");
    error = spg_get_error_code();
    printf("%s\n", spg_get_error_message(error));
  }
}
